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We complete our recently introduced theoretical framework treating the double quantum dot 
system with a generalized form of Hubbard model. The effects of all quantum parameters involved 
in our model on the charge stability diagram are discussed in detail. A general formulation of the 
microscopic theory is presented, and truncating at one orbital per site, we study the implication of 
different choices of the model confinement potential on the Hubbard parameters as well as the charge 
stability diagram. We calculate the charge stability diagram keeping three orbitals per site and find 
that the effect of additional higher-lying orbitals on the subspace with lowest-energy orbitals only 
can be regarded as a small renormalization of Hubbard parameters, thereby justifying our practice 
of keeping only the lowest-orbital in all other calculations. The role of the harmonic oscillator 
frequency in the implementation of the Gaussian model potential is discussed, and the effect of an 
external magnetic field is identified to be similar to choosing a more localized electron wave function 
in microscopic calculations. The full matrix form of the Hamiltonian including all possible exchange 
terms, and several peculiar charge stability diagrams due to unphysical parameters are presented in 
the appendix, thus emphasizing the critical importance of a reliable microscopic model in obtaining 
the system parameters defining the Hamiltonian. 
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I. INTRODUCTION 

The idea of quantum computation, introduced two 
decades agoP has attracted intense research interest be- 
cause of its ability to provide novel solutions to cer- 
tain problems that have been deemed intractable on a 
classical computerP The fundamental unit of a quan- 
tum computer is the qubit, which is implemented using 
quantum two-level (and sometimes multi-level) systems. 
Many physical systems have been put forward as candi- 
dates on which qubits can be embedded, such as nuclear 
magnetic resonance, 3 cavity QED, 4 and trapped ionsP 
Among these proposals, the solid state systems based on 
semiconductor quantum dots stand out as most promis- 
ing for the implementation of qubit since modern 
semiconductor industry allows great scalability for these 
systemsP^ffi2 Long coherence times have been experi- 
mentally achieved in both GaAtPEsl anc j gpHUl semi . 
conductor devices. One-qubit and two-qubit manipula- 
tions ha ve b een demonstrated in the GaAs quantum dot 
systems 12122103 i n silicon quantum dot systems, despite 
difficulties such as the fabrication-induced disorder^ and 
the valley de gree of freedom J22HH] high tunability has 
been reporte d 32 * 33 ! as a milestone toward coherent con- 
trol of the qubit. 

Two main effects are invariably present in all these 
experiments involving semiconductor quantum dots: the 
classical Coulomb interaction and the quantum fluctu- 
ations. On one hand, due to Coulomb blockade} 3 ^ the 
electron configurations in the quantum dot system are 
precisely controlled by electrostatic potentials (gate volt- 
ages). The details of electron configurations can be ex- 
tracted by an adjacent quantum point conta ct,^ and vi- 
sualized as the charge stability diag ram. 36 37 On the other 



hand, the quan tum effect, in particular the exchange 
interaction] 38 ! 39 ! plays an essential role in the qubit 
manipulation. The exchange interaction constitutes 
two-qubit operations in the Loss-DiVincenzo proposaP 
(which uses single-electron spin- up/down states as the 
qubit). In the double dot singlet/triplet proposal, the ex- 
change interaction controls the singlet-triplet level split- 
ting thus, along with an inhomogeneous magnetic field, 
provides arbitrary rotation along the Blo ch sphere and 
thereby achieves full single-qubit control P 13 * 23 ! 40 ! In a 
three-dot scheme the need for an inhomogeneous mag- 
netic field is eliminated, and the exchange interaction 
alone suffices for universal quantum computation. 

The charge stability diagram is fundamental for spin- 
based quantum computation: it is the starting point for 
all subs equen t qubit manipulation^^ and information 
readout)- 44 ! 45 ! Previous studies of the charge stability di- 
agram have focused on the classical Coulomb effects us- 
ing primarily the capacitance circuit modelJ^HSl which 
parametrizes the on-site and inter-site Coulomb inter- 
action as capacitances. The charge stability diagram is 
found through the electrostatic energy for a given elec- 
tron configuration on the double dot system. Although 
this classical theory has satisfactorily explained many 
aspects of experiments (especially the cases where the 
quantum fluctuations are weak), there are cases where 
the quantum effects, inevitably intertwining with the 
Coulomb interaction, lead to (sometimes substantial) de- 
formation of the charge stability diagra m from th at pre- 
dicted by the capacitance circuit modelF^ 33 1 50 1 51 ! A the- 
ory capable of reconciling both the classical and quan- 
tum effects becomes necessary, which will not only help 
us understand the quantum aspect of the charge stability 
diagram, but will also substantially extend its usefulness. 
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Attempts to explain the quantum aspects of the charge 
stability diagram have been many. Most of them have 
employed the quantum-mechanical two-level model to 
study the tunn eling of a single electron from one dot to 
the other EZHSHUl From the probability crossover of the 
two eigenstates (sometimes noted as the "excess charge" ) 
one can reliably extract the parameter for quantum- 
mechanical electron hopping (or tunnel coupling) t. A 
somewhat similar worlP2 has studied the influence of the 
microscopically calculated exchange interaction on the 
charge stability diagram. On the other hand, interest 
in applying Hubbard-like strongly correlated models to 
the quantum dot systenP^Ml arises in contexts related to 
the collective Coulomb blockade. In Ref. |59l a Hamil- 
tonian whose form is essentially very similar to our pro- 
posed generalized Hubbard model 50 has been derived us- 
ing the so-called "pocket-state" method. However, the 
power of Hubbard-like models to describe the charge sta- 
bility diagram h as lon g been neglected until Gaudreau 
and collaborator d 60 * 61 1 have applied the Hubbard model 
to triple, but not double, quantum dot system. 

In a previous publication,^ we have introduced a gen- 
eralized Hubbard model as a quantum generalization of 
the capacitance circuit model. We have shown that, 
with all quantum fluctuations suppressed, the generalized 
Hubbard model becomes the extended Hubbard model, 
which can be mapped to the classical capacitance circuit 
model exactly. The main advantage of our approach is 
that all terms allowed by symmetry arguments are kept 
in the Hamiltonian, which naturally accommodates all 
possible kinds of quantum fluctuations. We have further 
recognized that these quantum effects cannot be included 
in the capacitance model in any ad hoc manner. Be- 
cause the electron occupancies on each dot are no longer 
good quantum numbers, it breaks the basic assumption 
of the capacitance model that the electrostatic energy of 
the system is expressed in terms of electron occupancies 
on each dot, which are assumed to be integer numbers. 
Therefore, our generalized Hubbard model can be viewed 
as the quantum generalization of the classical capacitance 
circuit model, with the individual electron occupancy on 
each dot being a fluctuating quantum variable instead of 
being a fixed number as in the capacitance model. 

The quantum fluctuations may perturb the charge sta- 
bility diagram in a substantial way. In Ref. [50] we dis- 
cussed in particular the rounding effect of the boundary 
lines due to the electron hopping (tunnel coupling) t. A 
microscopic calculation was outlined to control the pa- 
rameters in the generalized Hubbard model, and we have 
shown how the charge stability diagram changes with 
the variation of quantum fluctuations induced by the de- 
formation of the microscopic confinement potential. In 
Ref. IHI] we applied o ur the ory to quantitatively explain 
two set of experiment d 32 1 33 1 on the silicon system and cal- 
culated, in particular, the tunnel coupling as a function 
of the height of the central potential barrier. The results 
are found to well describe the experimental data after 
appropriately choosing model parameters. 



In this paper, we address problems that have not been 
fully explored in our previous publications! 5 - ^ 5 -^ The pri- 
mary goal of this paper is to provide a comprehensive pic- 
ture of the theoretical framework that we have introduced 
in Refs. 1501511 thereby completing this series of study on 
the charge stability diagram of the double quantum dot 
system. In Refs. 1501511 we have proposed that the quan- 
tum fluctuations are essential for a in-depth understand- 
ing of the charge stability diagram, and have focused our 
attention on the tunnel coupling as it is manifestly the 
most significant quantum parameter. However, other pa- 
rameters such as the spin- exchange, pair-hopping and the 
occupation-modulated hopping are present as well, al- 
beit with smaller amplitudes. One of the central aims of 
this paper is to quantitatively estimate the effect of all 
these quantum parameters, which define the full quan- 
tum Hubbard model, on the charge stability diagram. 
Moreover, although we have outlined the microscopic the- 
ory and performed calculation in some specific cases in 
Refs. l5QI5"Tl there are several problems that have not been 
sufficiently covered in the previous publications. First, 
Ref. [5U] has focused on the biquadratic model potential 
and, since this model is rather special and involves only 
a few parameters, a dimensionless quantity r\ combining 
the height of central potential barrier and the energy of 
electron states is identified to quantify the quantum fluc- 
tuations and completely determine the geometry of the 
charge stability diagram. In Ref. EU the results calcu- 
lated from both the biquadratic model and the Gaussian 
model confinement potentials have been compared to ex- 
periments side-by-side. However a direct comparison of 
the two models in terms of the Hubbard parameters is 
missing. Although experimentally the electrostatic po- 
tentials dominate quantitatively, the detailed form of the 
actual confinement potential is obviously unknown. Thus 
it is important to understand the consequences of differ- 
ent choices of model potentials, under similar electro- 
static situations, on the Hubbard parameters as well as 
the charge stability diagram. We will therefore compare 
two models, biquadratic and Gaussian, in this paper with 
respect to the charge stability diagram. Seco nd, al l calcu- 
lations presented in our previous publications^^ as well 
as most of this paper are done assuming that each dot 
contains a single orbital and is only allowed to hold up to 
two electrons. We will discuss the case where this con- 
straint is lifted to allow three orbitals per site, and pro- 
jecting to the single-orbital subspace we study the effect 
of additional higher-lying orbitals on this single-orbital 
subspace. Last, we will also discuss the role of the mi- 
croscopic harmonic oscillator frequency in the Gaussian 
model potential, as well as the effect of external mag- 
netic fields. All of the se are new and not discussed in our 
previous publications! 50 * 51 ^ 

The remainder of this paper is organized as follows. In 
Sec. [Tl] we present a general microscopic formulation of 
the problem. In Sec. |III| we discuss in detail the general- 
ized Hubbard model proposed in Ref. [301 In Sec. |IV| we 
show charge stability diagrams calculated directly from 
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the generalized Hubbard model and discuss the effect of 
various parameters on the stability diagram. In Sec.[V|we 
present a very detailed discussion on the microscopic cal- 
culation, including the effects of different choices of con- 
finement potential, the additional higher-lying orbitals, 
the role of the harmonic oscillator frequency in the Gaus- 
sian model potential and the effect of the magnetic field. 
Sec. I VII is the conclusion. The full matrix form of the 
Hamiltonian and several cases with extreme (and possi- 
bly unrealistic) parameters are shown in the appendices. 

II. MICROSCOPIC THEORY: GENERAL 
FORMULATION 

We start with the general Hamiltonian of a system with 
N electrons, which consists of three parts: 

H{N) = H + H C +H Z . (1) 

Hq is the sum of the single-particle Hamiltonian for each 
electron, H n = J2iLi h {vi). He is the Coulomb interac- 
tion between electrons, 

JV-l N , 2 

i— i 

where k = 1/(47T£ £)- Hz is the Zeeman energy, 

N 

Hz=g*VBBY J S z l =E B S Z , (3) 

i=l 

where E B = g*(J- B B, S z = J2iLi S i- We neglect the 
spin-orbit couplings for simplicity. We further note that 
in our work, we do not include the interaction between 
the electrons on the quantum dots and the environment. 
Therefore Eq. ([!]) is understood as describing the elec- 
trons on the quantum dots only. In fact, many en- 
vironmental factors contribute to the decoherence pro- 
cess of the quantum dot system, such as the hyper fme 
interaction with the nuclear spin 

bath j 40 | 42 | 43 | 62 | 6 3] thg 

coupling to background impuritie d 64 * 65 ! and the phonon 
modes EMS] i n ^ n j g wor k ; we concentrate on the electronic 
interaction in the quantum dot system and highlight its 
consequence on the charge stability diagram. Therefore 
at this early stage of the theory we completely neglect 
the coupling to the environment. The elucidation of this 
problem is important in future studies. 

The single-particle Hamiltonian h (r) can be written 

as 

h(r) = ^[p-eA(r)] 2 + V(r) 1 (4) 

where m* is the effective mass of the electrons and V(r) 
is the potential (typically double-welled in the case of 
the triplet/singlet spin qubit) confining electrons on the 
xy-plane. We allow for a magnetic field B applied along 



the z axis, which couples to the electrons via the vector 
potential A. Here we have assumed that all electrons are 
experiencing the same magnetic field (in actual experi- 
ments this may not be the case). 

There are different choices of the detailed form of the 
confinement potential, but we assume that it is approxi- 
mately parabolic around the minima of the potential well 
so that the ground state single-particle wave functions 
are harmonic oscillator states. In this case, the single- 
particle Hamiltonian h (r) can be written as the sum 
of the Fock-Darwin Hamiltonian Hpd j and some non- 
harmonic potential Wj (r), h (r) = -Sfdj + Wj (r), with 

Hfdj = ^[ P -eA (r)] 2 + \m*uj) (r - v 3 f , (5) 

W i {r) = V{r)-\rn*d]{r-r i f . (6) 

Here, Uj is the corresponding harmonic oscillator fre- 
quency for dot j. We denote r'j = r — rj and express 
r'j as (r'j cos Op r'j sin 6 j) in Cartesian coordinates. The 
cigenfunctions of the Fock-Darwin Hamiltonian at the 
jth dot are the Fock-Darwin states satisfying 

H FVij \(pj, ni , mj (r'^O'j)) = e nj+inj _ \<p s , nj>mj (rj,0j)) . 

(7) 

where 

^n j+ , nj - = {n ]+ + 1/2) + huj- (rij- + 1/2) , (8) 
f],n 3 , mj (rpO'j) = exp [-i (eB/2H) (xjy - yjx)] 

■'■„, (V/ )/,V .,„ (r'j), (9) 

with 

i/j mj (6'j) = e im ^/V2K (10) 

and 

■«pW/^^'(|)- (id 

The index rij = 0, 1, • ■ ■ is the principal quantum num- 
ber, rrij = —rij, —rij + 2, ■ ■ • , rij — 2, rij is the azimuthal 
quantum number, and n T j = (rij — \rrij\) /2 is the ra- 
dial quantum number. The pairs of quantum numbers 
(rij, rrij) and (rij+,Tij-) are related by rij = rij- + %•+, 

rrij = rij- — rij + . denotes the Laguerre polynomi- 

als, l j0 = y/h/eB j + uj] /cj2 , and 

w,± = y/u] + o;2/4 ± w c /2, (12) 

where the Larmor frequency uj c — eB/m*. 

These Fock-Darwin states can be used either to solve 
the Hamiltonian approximately, or to generate an equiv- 
alent tight-binding model. In Ref. |SD1 we have obtained 
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a generalized Hubbard model using the Fock-Darwin 
states. There are infinitely many Fock-Darwin states, 
and it is impossible and unnecessary to keep all of them. 
For the M-dot system, we truncate the Fock-Darwin ba- 
sis in each dot to v levels, i.e., \<pi) = \^Pj jn ), where j labels 
quantum dots (1 < j < M) and n denotes the number of 
energy levels that we have kept (1 < n < v). We notice 
that the Fock-Darwin states within one dot are orthogo- 
nal to each other, (ipj, n \ Pj.n') = 0, but those from differ- 
ent dots are in general not orthogonal: (<Pj, n \ ¥>j,n') 7^ 0. 
We therefore build a new set of orthogonal bases by mak- 
ing the transformation^ 



(1*0 |* 2 ) •• 
= o- 1 '*Q< Pl ) 



\*Mv)Y 



\V>Mu)) 



(13) 



Here O is the overlap matrix (On* = (ipi\ fi')) generated 
by Fock-Darwin states in the single particle subspace. 
The new basis I*;) actually corresponds to, in a second- 
quantized form, c, |0) where ct creates an electron on 
site/orbital I with spin a. After the orthogonalization, 
the second quantized form of H and He can be built, 
in a standard way, as 



H'c 



X [ G h]hMM c h,t C h,l Cl 3^ Cl *^ 



^1 1^2 ^3 ^4 

.(2) 



The coupling parameters are 



t J 



o- C l 2 ,<T C h,<? c h,a 



F lul2 = J dr*l(r)h(r)V h (r) , 



(14) 



(15) 



(16) 



G 



fee 2 



ri - r 2 \ 



G 



(2) 

^1^2,^3)^4 



*i,(ri)* U (r 2 ), (17) 

fee 2 



-r 2 



(18) 



Then the full Hamiltonian can be written down in the 
second quantization form as H' = H' n + H' c + Hz ■ 

We note that the difference between our configuration 
interaction method and the traditional molecular orbital 
method-^! is that here we have transformed the eigen- 
value problem with non-orthogonal basis to a standard 
eigenvalue problem with orthogonal basis. 



III. GENERALIZED HUBBARD MODEL 

We now consider the specific case of a coupled double 
quantum dot system with each dot capable of holding 



up to two electrons. Effectively this means that we keep 
the s-orbital only, i.e. M = 2 and v = 1 in Sec. |Hj The 
case w here t he p-orbital comes into play will be discussed 
in Sec. VB As discussed in our previous publication,^ 
the system can be described by a generalized form of 
the Hubbard model which retains all terms allowed by 
symmetry: the total particle number N and the total spin 
S z are conserved. The one-body part of the Hamiltonian 
can be written as: 



H' 



AW 



tc\ a C2cr + H.C. 



(19) 



where = c? ia Ci a , fa denotes the chemical potential (or 
level energy) on site i (i — 1,2), and t denotes the inter- 
site hopping (or tunnel coupling). The Zeeman term is 

E B 

h z = -j- ( n n ~ n n + n n ~ n n) • ( 20 ) 

The two-body part H' c in general includes a Coulomb 
repulsion term H\j\ 

H LJ = Uin n n n + U 2 n 2t n 2i + U 12 (n n n 2i + n n n 2t ) 
+ (U12 - Je)(ni-fn 2t + n n n 2i ), 

(21) 

and a term Hj including spin-exchange (denoted by 
J e ), pair- hopping^2H25l (d enote d by J p ), and occupation- 
modulated hopping term^MISI (denoted by J t ): 



H, 



J ti n ia c \-g c 2o + H.C. 



(22) 



In the case with non-zero t, U\ and U 2 but all other 
parameters {Ui 2 , J e , J p , Jtx, Jt 2 ) being zero, one recov- 
ers the form of the usual Hubbard modeP^ (on a lattice 
with only two sites). The inter-site Coulomb repulsion 
U\ 2 was introduced in the study of the charge ordering in 
strongly correlated materials-^ and the model including 
the U\ 2 term is usually termed as the "extended Hub- 
bard model" PS In our previous work we have mapped 
the widely used capacitance model to the t = case of 
the extended Hubbard model and have argued that the 
extended Hubbard model is the minimal model that ex- 
plains the experiment.^ 

The spin-exchange and pair-hopping interactiorP^HZU 
have been studied in the context of orbital-selective Mott 
transition.^ For atomic orbitals simple relations exist be- 
tween the parameters (e.g. for <i-orbitals in free space, 
Ux = U 2 = U, J e = J p = J, U 12 = U - 2J and due 
to orthogonality t — J t — 0). However, in the quantum 
dot system these relations need not hold since the system 
usually has much lower symmetry than the free space or 
a lattice. The occupation-modulated hopping term J t 
has not been considered much in the literature except in 
certain aspects of superconductivity.^ 

The full Hamiltonian is a 16 x 16 matrix for the four- 
electron two-dot system (16 = 4 2 since each of the two 
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dots allows four possible quantum states). Because the 
total electron number N and total spin S z are conserved, 
it appears in a block-diagonal form. The details are pre- 
sented in Appendix [A] For a given (/i2,Mi) one finds the 
ground state by diagonalizing the Hamiltonian matrix. 
The charge stability diagram, plotted on a plane with 
M2 and Hi as axes, shows how the ground state changes 
as fi2 and /ii are varied. Experimentally, the chemical 
potentials of the dots are controlled by the gate voltages 
Vr and Vl and the charge stability diagram is plotted on 
the Vr-Vl plane. 

In the case of t = J e = J p — J t \ — Jti — 0, the 
ground state can be labeled as (n\, n 2 ) since (the elec- 
tron occupancy on dot i) is a good quantum number. 
As one or more of the t and J parameters becomes fi- 
nite, rii ceases to be conserved and the ground state is a 
linear combination of (711,712) states. The phase bound- 
aries between blocks with different N and S z are clearly 
defined as before since there is no mixing between them. 
Within an (N,S Z ) block we label the charge stability dia- 
gram using the (ni, 712) state that dominates the ground 
state. For example, t mixes the (1,0) and (0,1) states. 
For /ii > /12 the (1,0) state is the majority state, and 
vice versa. We then regard the line [i\ = M2 as the phase 
separator between two regimes in which (1,0) or (0,1) 
are the majority states. Experimentally, these boundary 
lines sepa rating different states within the same N block 
are diffuse^EH due to the hybridization of those states. 
To fully reproduce this fading effect, a conductance cal- 
culation involving excited states is needed.^ In our work 
we concentrate on the ground state feature so the fading 
of the boundary lines is beyond the scope of this work. 

In this paper we assume Ux = U 2 = U (except in 
the discussion of Fig. [2]), and J t \ = J t2 = Jt- This 
means that there is a 1, 2-permutation symmetry and 
the charge stability diagram is symmetric with respect 
to the line Ml — M2- Under this assumption plus a con- 
dition that Jt = 0, the charge stability diagram has a 
particle-hole symmetry (symmetric with respect to the 
line Mi + fi2 = U + 2fi2 — J e ), as can be seen from 
the matrix form of the Hamiltonian in Appendix [X] In- 
clusion of Jt effectively changes the values of t in high- 
occupancy states, thus breaking the particle-hole symme- 
try. Of course allowing U\ 7^ U2 would lead the stability 
diagram not being symmetric with respect to the line 
Mi = A*2- 

The calculated stability diagram from the generalized 
Hubbard model will be discussed in the next section. 
Analysis of the stability diagram shows that some of the 
parameters can be extracted directly from experimen- 
tal plots. Alternatively, all the parameters of the model 



assumed to be linear combinations of (Vr, Vl)! 60 * 61 ^ 



can be calculated from the microscopic theory [Eqs. ( 16 1, 



(171, and ( 18 1] using a microscopic confinement potential 



model. Therefore, the generalized Hubbard model es- 
tablishes a quantitative correspondence between the mi- 
croscopic theory and the experiments. Since the experi- 
ments take Vr and Vl as the basic variable, it is useful 
to relate Vr, Vl to ji 2 and Mi- In the literature, (m2,Mi) 



Mi = |e|(aiVr, + /3iVR)+7i, 
M2 = |e| (a 2 Vr + /3 2 Vl) + 72; 



(23) 
(24) 



where 71 and 72 are constant energy shifts. 

In our previous worlP^ we have presented a map- 
ping between the capacitance model and the extended 
Hubbard model, which suggests that the coefficients in 
Eqs. (23 1 and (24) are related to the parameters bjE^ 



«i 



(U 2 - U 12 ) Ui 

U!U 2 -U? 2 ' 



OL2 



(Ui ~ U12) U2 
UxU 2 - UI2 



(25) 



and /3i.2 = 1 — ai.2- When the two dots are symmetric, 
we define a = a\^ 2 , ft = flx 2 and Eq. |25|) reduces tcP^ 



U 



U + Ut 



(26) 



and /3 = 1 — a. 

We also note that the relation implied by Eqs. (23) 
and ( 24 ) is not a unitary transformation. Since we are 



primarily interested in the lengths of segments on lines 
fJ-i—fJ-2 = const., the lengths do conserve, up to a factor of 
electron charge e, upon transforming (/i 2 , Mi) ~ ^ (Vr, Vl) 
in our theory. However, one must be cautious in fitting 
the experiments^ since the experimental values of the 
gate voltages should sometimes be rescaled by a factor 
in order to apply our theory. This uncertainty between 
the parameter sets (^2,^1) and (Vr, Vl) is unavoidable 
within the scope of our theory and can only be resolved 
with precise experimental information. 



IV. RESULTS FROM THE GENERALIZED 
HUBBARD MODEL 

In this section we present charge stability diagrams 
calculated directly from the generalized Hubbard model. 
The parameters U, U\ 2 , t, J e , J p , and J t are assumed 
to be independent of fj, 2 and fix- This is different from 
the calculations from the microscopic theory where the 
overlap integrals of wave functions at different locations 
on the stability diagram are not guaranteed to be the 
same and the parameters indeed change slightly at dif- 
ferent locations on the charge stability diagram. We also 
assume Eg — in this section for simplicity. Throughout 
this section we discuss the effect of various parameters in 
the generalized Hubbard model on the charge stability 
diagram. To facilitate the comparison we plot all figures 
within the same /12, Mi range. We also note that although 
we set the unit of parameters to be meV (so that the re- 
sults are comparable with experiments and the results 
from microscopic theory), it can essentially be arbitrar- 
ily chosen with the actual energetics being determined 
by the details of microscopic confinement. 

The simplest case is the system of two decoupled quan- 
tum dots with on-site Coulomb interaction only, i.e. 
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FIG. 1: Charge stability diagram calculated at U = 6.1 meV, 
U12 = 2.5 meV, and t = J e = J p — J t = 0. 



FIG. 2: Charge stability diagram calculated at Ui = 5 meV, 
f/ 2 = 7 meV, U 12 = 2.5 meV, and t = J e = J p = J t = 0. 



U12 = t = J e = Jp = J t = 0. Its charge stability di- 
agram has a checkerboard pattern which obviously does 
not fit the experiment. As mentioned in Ref. 1501 the ex- 
tended Hubbard model (with hopping terms neglected) is 
the minimal model that explains the experiment. There- 
fore we start with the case < U\2 < U, while keeping 
t = J e = J p = Jt = 0. A typical calculation is shown in 
Fig. [T] The parameters U and U\2 are chosen to be the 
same as Fig. 1 of Ref. [50] and the plot is identical to the 
solid line (i = result) of Fig. 1(b) in Ref. [SpJ 

We examine the charge stability diagram in detail, in 
particular we are interested in the form of the phase 
boundaries and the coordinates of the "triple points" p3 
defined as the points on the diagram neighboring three 
different phases. In this case, all eigenenergies are lin- 
early dependent on jUj, U, and U\ 2 . Therefore all phase 
boundaries are straight lines, which can be seen clearly 
from Fig. [T] The triple points are lettered on the fig- 
ure with notations that exhibit symmetry considera- 
tions. The 1,2-permutation symmetry does two things: 
First it ensures that triple-points A, B, A', B' lie on 
the line \i\ — fi 2 . Second it implies that the coordi- 
nates of C and C can be found by interchanging indi- 
cies 1 and 2 in the corresponding coordinates of C and 
C . The particle-hole symmetry implies that fi^iX') = 
U + 2U 12 - J e -fJ,2,i(X), where X = A,B,C,C. 

The precise coordinates (/X2,/xi) of the triple points 
can be readily calculated. We have A(0, 0), B(Ui 2 , U12), 
C(U, U12), and C'(U + U 12 ,2U 12 ). (Coordinates of other 
points can be found by symmetry.) We notice that the 
length of the line segments AB = CC = \/2Ui 2 and 
BC = U. In particular, the length of AB and CC does 
not change upon transforming from (/Lt2,/iti) to (Vr, Vj,). 
Therefore from the experimentally measured charge sta- 
bility diagramj22l one can read off the value of U12 from 
the length of the phase boundary between (1, 1) and 
(0,2), as what we have done in Ref. [50] In Ref. l50l we 
have also extracted the value of U from the slope of the 
phase boundary between (1,1) and (1,2), according to 



Eq. (26 1. We note that the readout of U12 is implied 



by the extended Hubbard model per se, while the value 
of U requires knowledge of the relation between (^2, /ii) 
and (Vr, Vl). While the precise form of this relation is 
in general unknown, the mapping between the capaci- 
tance model and the Hubbard model provides one in the 
simple cases where all quantum fluctuations vanish. In 
the cases where the mapping is not necessarily valid, one 
must exercise caution. 

The case < U\2 < U considered here is physically 
reasonable, which can be confirmed by calculations using 
the microscopic theory (see, e.g. Fig. [9]). For complete- 
ness we also consider the case with U12 > U and the 
results are shown in Appendix [Bj 

In Fig. [2] we show the result calculated allowing 
U\ 7^ U 2 . Both 1,2-permutation symmetry and particle- 
hole symmetry are broken. We have the coordinates 
(/^2,/ii) of the triple points as A(0, 0), B(Ui2, U12), 
C(U 2 , U 12 ), C'(U 2 +U 12 , 2U 12 ), C{U 12 , Ui), C 7 (2U 12> U 1 + 
U 12 ), B'(U 2 +U 12 , Ux+Uia), and A'(U 2 +2U_ 12 , Ui+2U 12 ). 
Note that the lengths of line segments AB and CC re- 
main v / 2t r i2- Since the experimental plot of Ref. |2"51 does 
not show the whole range of the stability diagram, the 
explanation of U± 7^ U2 is still consistent with the ex- 
periment, i.e., an alternative fit allowing U± ^ U2 is still 
valid. 

Fig. [3] presents a typical result calculated with a finite 
t while J e = J p = J t = 0. The data shown here is iden- 
tical to that shown as the red dotted line in Fig. 1(b) of 
Ref. [SOI The main effect of t is to mix states with the 
same total electron number N, except N = 0, 4. In this 
case one needs to identify the dominant (711,712) state in 
the true ground state to label the corresponding region 
on the charge stability diagram, as mentioned above. In 
Fig. [3] we use a star as a superscript of (ni,n 2 ) to in- 
dicate that this is the majority state in a linear com- 
bination rather than the true ground state. Most of the 
phase boundaries are curved and the one separating (0, 0) 
and (1,0)/(0,1) complex can be identified as hyperbola 
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FIG. 3: Charge stability diagram calculated at t — 0.6 meV, 
U = 6.1 meV, U 12 = 2.5 meV, J e = J p = J t = 0. 



= i 2 - [The boundary between (2, 2) and (2, 1)/(1, 2) 
complex can be found by particle-hole symmetry condi- 
tion.] The singlet combination of the (1, 1) states domi- 
nates the (1,1) component of the stability diagram [de- 
noted by (1, l)!?* in the figure], whose energy eigenvalue 
can be found by diagonalizing Eq. ( A-6 ) where simple 
analytical formula does not exist. The coordinates of A 
and B can be found analytically: [12(A) = Hi(A) = —t, 
fi 2 {B) = miB) = i+(C/+C/ 12 )/2- v /4i 2 + (U - (7 12 ) 2 /4, 
while the coordinates of C and C cannot be found ana- 
lytically. It is evident, however, that CC is stretched by 
the introduction of tunnel coupling t. 

To understand the effect of J e we first switch off t. 
Fig. pja) shows a typical charge stability diagram calcu- 
lated at < J e < U\2 with t = J p = J t — 0. Since 
the triplet combination of (1,1) state [labeled by (1, 1)T] 
decouples from Eq. (A-5), all eigenenergies are again lin- 
ear functions of U, U12, and J e thus the phase bound- 
aries are linear. The coordinates (^2,^1) °f t ne triple 
points are A(0, 0), B(U 12 - J e , U 12 - J e ), C(U, U 12 - J e ), 
and C(U + U12, W\2 — Je)- The length of line segments 
~AB = V2(U 12 - J e ), and CC 7 = V2U U . If the whole 
charge stability diagram is measured, one can read off the 
value of J e immediately from the difference between the 



two segments AB and CC (the effect of other parameters 
t, J t must be assumed to be small, though). Typically J e 
is smaller than U or U12 by at least an order of magni- 
tude. Therefore the J e = 1 meV (with U = 6.1 meV, 
U12 = 2.5 meV) shown in Fig. Qa) is slightly exagger- 
ated in order to make its effect clear. In actual systems 
one expects the effect of J e on the charge stability dia- 
gram to be very small. For a brief discussion of the case 
Je > U12, see Appendix [SJ 

The pair hopping interaction mixes the (2, 0) and (0, 2) 
states. If t + Jt 7^ 0, it in addition mixes the (1, 1) sin- 
glet with the (2,0) and (0,2) states, as can be seen from 
Eq. ( TO} . Fig.gb) shows a case with a finite but small 
J p . Comparison to Fig.Qa) reveals that the change im- 
posed by J p is very small. The boundaries of (0,2)* 
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FIG. 4: Charge stability diagram calculated at J e = 1 meV 
and (a) J p = 0, (b) J p = 1 meV. U = 6.1 meV, U 12 = 2.5 
meV. t = J t = 0. 
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FIG. 5: Charge stability diagram calculated at U = 6.1 meV, 
U12 = 2.5 meV, t = 0.6 meV, J e = J p = 0.2 meV. J t = 0. 



and (2,0)* region are indeed curved but the curvature 
is hardly detectable. The boundaries between N = 0, 1 
blocks and N — 3, 4 blocks remain the same. The cor- 
rection to the coordinates of the triple points C and C 
(and their mirror C and C) can be shown to be of order 
Jl/{U-U l2 ). 
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FIG. 6: Charge stability diagram calculated at (a) J t — (b) 
Jt = 0.5 meV and U = 6.1 meV, U w = 2.5 meV, t = 0.6 
meV, J e = J p = 0.1 meV. 



Comparison of Fig. [3] and Fig. [4^a,) shows that the 
(1,1) component of the stability diagram has interesting 
behavior: for t > J e the singlet combination dominates 
the (1,1) block, while for t <C J e the triplet combination 
domin ates. This can be understood from the eigenvalues 
of Eq. ( A-5 ). The triplet has energy — fii — fi 2 + U\% — J e 
while the singlet has energy —fix — /12 + U12 + J e phis 
contributions due to t + J t . When t + J t is small the 



triplet has an energy 2 J e lower than the singlet, thus be- 
coming the ground state. However when t + J t is large 
the correction due to t + Jt exceeds 2J e , making the sin- 
glet the dominating ground state. One expects an inter- 
mediate value of (t + J t )/J e such that the singlet and 
triplet co-exists. This situation is shown in Fig. [5] where 
the boundaries separating the singlet and the triplet are 
plotted as dashed lines. The condition that the singlet 
state dominates can be analytically obtained as: 



1 



(t + JtY > MJe + J P + U- U 12 ), 



(27) 



which is typically satisfied in most microscopic calcula- 
tions of double quantum dots. The condition that the 
triplet dominates without having singlet regimes on its 
sides cannot be expressed analytically. 



At the end of this section we turn on all parameters 
that we have discussed and consider the effect of J t . 
Fig. |6|a) shows a case with all parameters being non- 
zero except Jt- The figure is essentially very similar to 
Fig. [3] since the values of J e and J p are very small and 
have Tittle effect. The particle-hole symmetry is evident. 
Fig. ^b) shows the effect of a relatively large J t - First, 
the N = 0, 1 parts of the stability diagram remain un- 
changed. Second, the N = 2 parts [(0, 2)/(l, 1)5/(2, 0) 
complex] are slightly stretched towards the N = 3 block. 
The N = 3 block [(2,1)/(1,2) complex] is deformed 
and both of its boundaries are substantially rounded and 
smoothened. We expect that this is due to the change 
of an effective t to t + J t (N = 2) or t + 2J t (N = 3). 
Overall, it is clear that J t breaks the particle-hole sym- 
metry, making the boundaries between high-occupancy 
states smoother. 



V. RESULTS FROM THE MICROSCOPIC 
THEORY 

As disc ussed i n Refs. \EU\ and [STJ the microscopic 
calculatiorPsEHUHlU is required to constrain the param- 
eters of the Hubbard model in the physically relevant 
regime. In other words, although the Hubbard model 
by itself can have arbitrary parameters, a given physical 
double-dot system is restricted by the realistic confine- 
ment potential which would severely restrict the physi- 
cal Hubbard parameters for the system. The details of 
the application of microscopic theory to our model are 
described in Sec. [H] In this section we discuss the ap- 
plication of the microscopic theory to our problem. We 
shall primarily use two different models of potential: a bi- 



quadratic form [Eq. ( 28 1] and a Gaussian form [Eq. (29)]. 
Both of these model potentials are reasonably realistic 
and used extensively in describing d ouble-do t systems. 
The biquadratic potential is defined aipSMEEl 



(x + a) 



V Q (x,y) = Mm{ n fl[(. 

^[(,- a) 2 + y 2 ]-M 2 }, 



Mi) 



(28) 



(note that here the two dots are assumed to be symmet- 
ric) and the Gaussian potential reads 



V G (x,y) = -Vi exp 



(x + a) 2 + y 2 



f 2 



-V2 exp 



jx-a) +y 2 
l 2 

'c/2 



Vb exp 



x 2 + y 2 " 1 



I 2 



(29) 



Note that Eqs. (28) and (29) have slightly different forms 



than that in Ret. 1511 in order to facilitate the comparison 
of the effect of different potentials on the charge stability 
diagram in Sec. V A This will be explained in Sec. |V A| 
Besides, we a lso c onsider the effects of the additional p- 
orbitals (Sec. VB), and the role of the harmonic oscilla- 



tor frequency in the calculation with the Gaussian model 
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FIG. 7: (Color online) Profiles of biquadratic (black solid 
line) and Gaussian (red dashed line) model potentials under 
two different values of detuning energy e = y,2 — fii. For the 
biquadratic potential, huio — 3.956 meV, a = 33 nm. The 
parameters of the Gaussian potential are derived from the 
biquadratic one. The gray dotted lines denote the energy 
reference values of the chemical potential. 



potential (part of Sec. VC). All calculations mentioned 



above are done without magnetic field, and we consider 
the influence of the magnetic field in Sec. V C 

The dielectric constant for GaAs e = 13.1 and the 
effective mass (m* = 0.067m e ) are used throughout the 
microscopic calculation. 



A. Influence of model confinement potentials 

In the microscopic calculation, the confinement poten- 
tials in the quantum dot systems are usually represented 
by models. Since one has several choices of models, it is 
useful to understand the consequences of different model 
potentials on the charge stability diagram. A similar 
comparison has been done for the exchange coupling in 
Ref. 1201 but in light of the generalized Hubbard model we 
would like to understand the influence on the overall ge- 
ometry of the charge stability diagram and the Hubbard 
parameters. We match the two models of potentials ac- 
cording to their minima and the central barrier (see Fig 
for examples at two different values of detuning energ 
because these quantities are directly related to the elec- 
trostatic potentials that one can experimentally control. 
The biquadratic model potential is chosen first because it 
contains only two parameters. The chemical potentials 
fj>i t 2 enter the problem in two ways: first, they appear 
in the formal definition of the model [Eq. (28)]; second, 
the wave function solution to the Schrodinger equation 
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FIG. 8: (Color online) Charge stability diagrams calculated 
from biquadratic [panels (a) and (b)] and Gaussian [panels 
(c) and (d)] model potential different inter-dot distances. The 
parameters for the biquadratic model are huo = 3.956 meV, 
a — 33 nm (panels (a) and (c)), and a = 28 nm (panels (b) 
and (d)). The white crosses [indicating the points at which 
the comparison of Hubbard parameters (Fig. [9J is conducted] 
are located at fii = fj,2 = 9.743 meV for panels (a) and (c), 
and /ii = jj.2 = 10.095 meV for panels (b) and (d). 
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FIG. 9: (Color online) Parameters in the generalized Hub- 
bard model as functions of half of the inter-dot distance a, 
calculated at the white crosses shown in Fig. [8] using both 
biquadratic and Gaussian potentials. 
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and ( 16 )] . In our Hund-Mulliken calculation for the bi- 



quadratic potential, the chemical potentials defined in 
these two ways are approximately equal; while for the 
Gaussian model the chemical potential values calculated 
from the wave functions are different from V± and V2 de- 
noted in Eq. ( 29 1 . Since we want to directly compare 



the results from the two models under the same electro- 
static situations, we tune parameters (Vi.,2) such that the 
calculated chemical potential values match. This is the 



main reason that Eq. ( 29 ) is written in a slightly differ- 



ent way than Ref. 1511 Since the Gaussian model contains 
more parameters than the biquadratic model, some of the 
parameters in the Gaussian model have to be fixed: we 
fix lb = 20nm. The reference energy values (shown as 
light gray dotted lines in Fig. [7| correspond to the zero 
of the chemical potential in Fig. [HJ A close examination 
of Fig. [7] reveals that the main difference between the two 
models is that the Gaussian model has narrower potential 
wells. This implies stronger on-site Coulomb interaction 
and weaker inter-dot tunnelling, as shall be seen in the 
following discussion. 

In Fig. [8] we compare the microscopically calculated 
charge stability diagram for different models of poten- 
tials and inter-dot distances. Note that all charge stabil- 
ity diagrams shown in this section have an overall energy 
shift relative to those shown in Sec. IIVI due to the zero 
point energy of the harmonic oscillator states. Panels (a) 
and (b) show the results calculated from the biquadratic 
potential, while the results from the Gaussian model are 
shown in panels (c) and (d). The value of a for panels 
(a) and (c) is a = 33 nm, while that for panels (b) and 
(d) is a = 28 nm. The parameters a and ujq (given in the 
caption) are for biquadratic model, while those for the 
Gaussian model are derived hereby as explained above. 
We see that the overall shapes are similar for all panels. 
Panels (c) and (d) are shown with a larger x- and y-range 
than panels (a) and (b), because the Gaussian model has 
a larger on-site Coulomb interaction than the biquadratic 
model at the same interdot distance, shifting the phase 
boundaries to higher energies. This is consistent with 
the qualitative argument above that the narrower poten- 
tial wells of the Gaussian model lead to larger on-site 
Coulomb interactions. Moreover, the rounding effects 
near the triple points in panels (b) and (d) are more pro- 
nounced than panels (a) and (c) , which originates from a 
smaller inter-dot distance 2a, leading to stronger quan- 
tum fluctuations. These arguments will also be confirmed 
in the calculated Hubbard parameters below. We have 
noted that for the biquadratic model, the chemical po- 
tential values derived from Eqs. ( |l~4"| and (16 1 are some- 
times slightly different from the values in Eq. |28~1). This 



implies that in Fig. 8ja) and (b) some of the boundary 
lines, which are parallel to the x- or j/-axis in the Hubbard 
model calculation, are slightly tilted relative to the axes. 
This effect is more pronounced for a smaller inter-dot dis- 
tance 2a [Fig. [8^b)], since a stronger overlap between the 
wave functions in the two dots leads to a larger devia- 
tion of the calculated chemical potential values from the 



values directly implemented in the model. This tilting 
effect is not present in the Gaussian model calculation 
[Fig. f8fc) and (d)], as the chemical potential values are 
chosen to be the ones calculated from the wave functions. 

We would also like to understand the Hubbard pa- 
rameters as functions of half the inter-dot distance a 
for the two models. We choose one particular point on 
each charge stability diagram in a consistent way to com- 
pare the parameters. These points are represented as the 
white crosses in the (1,1) components of the charge sta- 
bility diagrams in Fig. [SJ they are defined as the mid- way 
point of the (1,1) segment of line /ii = /i 2 for biquadratic 
model calculations with different values of a, and those 
coordinates are inherited in the Gaussian model calcula- 
tions, in spite of the fact that they are not in the center of 
(1,1) regime any more (which is, again, a consequence of 
the more pronounced localization effect of the Gaussian 
potential). In Fig. [8] the white crosses are actually at the 
same location for panels (a) and (c) , as for panels (b) and 
(d). We have calculated the Hubbard parameters at these 
points for a range of a values (which are bracketed by 
those shown in Fig. [8]) , and the results are summarized in 
Fig. [9j Fig. |9|a) reveals that the on-site Coulomb repul- 
sion U is independent of the inter-dot distance, indicating 
that the wave functions are well localized. U is larger in 
the Gaussian model than in the biquadratic model by 
about 35%, which is a consequence of narrower potential 
wells in the Gaussian model: the electron wave function 
is more localized than the biquadratic model. For the 
same reason, all parameters characterizing the inter-dot 
interactions are smaller in the Gaussian model than in the 
biquadratic model, as quantitatively shown in Fig. (9|b)- 
(f). Fig. [9^b) shows that the inter-dot Coulomb inter- 
action is only slightly decreased in the Gaussian model, 
while Fig. [9^c) and (e) show that the hopping (tunnel 
coupling) as well as the occupation-modulated hopping 
for the Gaussian model is decreased by more than 40% 
from that for the biquadratic model. This implies that 
the difference in the dimcnsionlcss ratio t/U can only be 
more pronounced, as shown in Fig. [9]f f ) . In Fig. [9^d), 
one sees that the spin-exchange and pair-hopping inter- 
action is substantially smaller in the Gaussian model: it 
is decreased by an order of magnitude from that for the 
biquadratic model. 

Our results indicate that in modeling the double quan- 
tum dot systems, the particular choice of model confine- 
ment potential leads to quantitative changes in the Hub- 
bard parameters, although the qualitative behaviors are 
similar. In a realistic study of these systems, the first- 
principles calculation of the exact form of the confine- 
ment potential from Poisson equation would be required. 
However, this demands higher precision in both experi- 
ments and theoretical calculations, which should be pos- 
sible in future studies if more precise information about 
the lithographic gates creating the confinement potential 
is available. 
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FIG. 10: (Color online) Charge stability diagrams calculated 
from the Gaussian potential using (a) s-orbitals only and (b) 
s- and p-orbitals. The parameters corresponding to the bi- 
quadratic potential are a — 30 nm and Hcjo = 3.956 meV. 
Panel (c) shows the full charge stability diagram including all 
49 components, while panel (b) shows the stability diagram 
in the two-electron regime. Panel (b) is a replica of the area 
enclosed by the white frame in the lower left corner of panel 
(c). 



B. Effect of additional p-orbitals 

All of the calculations shown in this paper so far, as 
well as those shown in Refs. |SD] and [5TJ bear the same 
feature that each quantum dot is only allowed to hold up 
to two electrons. This means that we have taken M = 2 
and v = 1 in the general discussion in Sec. [TTJ, keeping 
the s-orbital (denoting the lowest-lying state) only. In 
the discussion of silicon devicesp^ a multi-electron multi- 
band situation is encountered. In Ref.EUwe have reduced 
the full multi-band problem to an effective two-electron 
regime and have argued that the multi-band effect can be 
regarded as a renormalization of Hubbard parameters. In 
this subsection, we examine this idea in a slightly differ- 
ent context: we consider the additional p-orbitals (which 
are doubly degenerate and can be labeled as p x and p y 
orbitals), the lowest excited states except the s-orbitals 
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FIG. 11: (Color online) Effective on-site Coulomb interaction 
[panel (a)] and inter-site Coulomb interaction [panel (b)] of 
the results shown in Fig. [10] as a function of the detuning 
energy e = fj.2 — Hi calculated keeping s-orbitals only (black 
solid line) and s- and p- orbitals (red dashed line). The traces 
where calculations take place are shown as white diagonal 
lines in Fig. 10 a) and (b). As s is swept from the top left 
corner to the bottom right corner, the dominating ground 
state is changed from (2,0) to (1,1), and then to (0,2). The 
black dotted lines and red dash-dotted lines denote the phase 
boundaries (where the probabilities of states cross) in the case 
of s-orbitals only and s- and p-orbitals, respectively. 



considered. This can be regarded as taking M — 2 but 
v = 3 in Sec. [TT] The full charge stability diagram has 
49 components ranging from (0,0) to (6,6). We focus 
on the part with population up to two electrons per dot 
(2,2), and study how some of the Hubbard parameters 
are renormalized by the additional orbitals. 

Fig. [10] compares the charge stability diagram calcu- 
lated keeping s-orbitals only and that calculated keeping 
both s- and p-orbitals. Panel (a) is the result of keeping 
the s-orbitals only. Panel (c) shows the full charge sta- 
bility diagram for the case with both s- and p-orbitals, 
including all 49 components. Its lower left corner corre- 
sponds to the two-electron regime and we replicate the 
area enclosed by the white frame as panel (b) and com- 
pare that to panel (a) side-by-side. Strong similarity be- 
tween panel (a) and panel (b) is evident and the differ- 
ences are only in the details. In Fig. 10 'b), the stability 



diagram are slightly shifted toward the origin, and the 
diamond of the (1,1) component shrinks its size. This 
means that the Coulomb interaction U is renormalized 
to a smaller value by the additional orbitals. 

To quantify this idea we study the "effective" on-site 
Coulomb interaction (U cS ) and inter-site Coulomb inter- 



in Fig. 
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action (C/ff ) 

effective Coulomb interactions is that we want to intro- 



The purpose of studying the 
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duce a method which may be useful for a reduction of the 
general multi-band problem to the two-electron regime. 
For the case studied here (retaining three orbitals) it is 
still possible to calculate all parameters in the general- 
ized Hubbard model, but for a general multi-band prob- 
lem this calculation is cumbersome and unnecessary, as 
all the qubit manipulations are done in the two-electron 
regime and it is sufficient to study the effect of the addi- 
tional band on this regime. 

For a general multi-band problem, the total electron 
number N and total spin S z should still be conserved. 
The full Hamiltonian is block-diagonalized as per {TV, S z } 
and we denote the n th lowest energy eigenvalues in each 
block as Ei N ' S ' } . Then when the (2,0) or (0,2) is the 
dominant ground state, the effective on-site Coulomb in- 
teraction can be defined as the energy difference between 
the doubly occupied state and the singly occupied state: 



XT 



off 
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2E 



{b±i} 



(30) 



Note that here we have used the fact that the two dots are 
symmetric. When they are asymmetric, Eq. (30) gives 



the Coulomb interaction on the first dot (U\) when cal- 
culated at chemical potential values such that (2,0) dom- 
inates the ground state, and it gives Ui when calculated 
in cases where the (0,2) dominates the ground state. 

Moreover, for cases where (1, 1) dominates the ground 
state, the effective inter-site Coulomb interaction C/°f is 
defined as the difference between the energy of the state 
with each dot holding one electron and the total energy 
of states that one electron is allowed to fill one and the 
other dot respectively. This means that in the {N — 
1,S Z = ±i} block both the lowest and the second lowest 
energy eigenvalues are involved, as we want the electron 
to occupy a different dot after one of the two dots is 
occupied. Then {Tff is expressed as 



cff 
2 



E {2 ' ±1} 



E 



{b±|} 



E, 



{i.±|} 



(31) 



In Fig.|lO[a) and (b) we have drawn two diagonal lines 
in the N = 2 block of the charge stability diagram, as 
this is the regime which attracts the most interest. Along 
these lines we have calculated the effective on-site and 
inter-site Coulomb interaction as functions of the detun- 
ing energy e = fi2 — Mi) an( i the results are shown in 
Fig. [TTJ The black solid lines in Fig. [TT] show the case 
with s-orbitals only, corresponding to Fig. [To]a), while 
the red dashed lines show the case with s- and p-orbitals, 
corresponding to Fig. [To](b) . One must be cautious when 
interpreting Fig. [TTJ as the definitions of the effective 
Coulomb interaction Eqs. p0| and (31) is only valid 



the 



for a selective range of parameters that some particu- 
lar states dominates the ground state. In Fig. 11 a) 
values shown in the range |e| > 10 mcV gives the effective 
on-site Coulomb interaction: For the case with s-orbital 
only, U cS — 10.38 meV. This is consistent with the value 
in the Hubbard model U = 10.38 meV. For the case with 
both s- and p-orbitals, U cS changes slightly with the de- 
tuning energy but is bound to a narrow range between 



9.24 and 9.68 meV. This indicates that the main effect of 
the p-orbital is a renormalization of the effective on-site 
Coulomb interaction to a slightly smaller value, which 
has a very weak dependence on chemical potential val- 
ues. The results shown in the range |e| < 10 meV are 
meaningless and will not be discussed here. Fig. [TT|b) 
shows the calculated an d the range \e\ < 6 meV is 
directly relevant here. For the case with s-orbital only, 
U^2 = 1-88 meV, consistent with the value in the Hub- 
bard model U\2 = 1.88 meV. For the case with both s- 
andp-orbitals, C/°f = 2.00 mcV. Therefore the additional 
p-orbitals does not change effective value of U\2 apprecia- 
bly. The calculated values of U eS and C/ff are also con- 
sistent with that extracted directly from Fig. [lO^a) and 
(b) . The main result in this subsection is that the effect of 
additional orbitals can be regarded as a small renormal- 
ization of the Hubbard parameters and the main physics 
is indeed captured in the two-site one-orbital Hubbard 
model. 

We note, that in Fig. [TT|b) the U 12 value extracted 
from the range \e\ > 12 meV is actually related to the 
Coulomb interaction between s and p orbitals (which can 
be denoted as U sp ), and this information is in principle 
encoded in the stability diagram of Fig. 10 [c) . In fact, the 
full characterization of the three-orbital charge stability 
diagram in Fig. |10[ c) shows rich behavior: the compo- 
nents of the stability diagram generically have smaller 
area for highly occupied states than the one-orbital sub- 
space we are considering, for example the (3,3) and 
(5,5) states spans a substantially smaller range on the 
charge stability diagram. This indicates that the on-site 
Coulomb interaction is weaker for higher-lying orbitals 
than the s orbitals. In fact, as the p orbitals are popu- 
lated, the effective central potential barrier between the 
two dots are effectively lower since p orbitals themselves 
have higher energies, and the corresponding wave func- 
tions are less localized. This naturally leads to the fact 
that on-site Coulomb interaction for higher-lying orbitals 
is relatively smaller. Moreover, through the experimen- 
tal observation of the charge stability diagram, one can 
in principle infer the effective microscopic confinement 
potential. The full understanding of this multi-orbital 
problem is, however, beyond the scope of the present 
work and is an interesting topic for future investigation. 



C. Harmonic oscillator frequency in the Gaussian 
model and the effect of external magnetic field 

In this subsection we study the role of the harmonic 
oscillator frequency in the Gaussian model and the effect 
of external magnetic field. 

There are intrinsic difficulties when applying the Gaus- 
sian model potential in the Hund-Mulliken approxima- 
tion: as the chemical potentials change their values, the 
two potential wells are deformed such that a harmonic- 
oscillator-type expansion of the confinement minima 
yields considerably different values of the harmonic os- 
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FIG. 12: (Color online) Charge stability diagrams calculated 
from Gaussian potential with a — 30 nm and huo — 3.956 
meV. (a) B — 0, the zero-point energies of Fock-Darwin states 
are chosen according to the microscopic confinement potential 
at /ii = fi2 = 9.938 meV (indicated by the white cross) with 
fu^oa = 6.487 meV. (b) B = 0, the zero-point energies of 
Fock-Darwin states are ftojoc = 9.731 meV. (c) B = 0, the 
harmonic oscillator frequencies of the Fock-Darwin states are 
chosen according to the microscopic potentials for each given 
Hi and jj.2- (d) B — 6 T, the zero-point energies of Fock- 
Darwin states are Hljog = 6.487 meV. 



cillator frequencies. Since the harmonic oscillator fre- 
quency is the central variable for the Fock-Darwin states, 
its variance leads to a change in the Hubbard parameters 
as well as deformation of the charge stability diagram. 
Moreover, the values of the chemical potentials them- 
selves rely on the calculated Fock-Darwin states, which 
further depend on how the harmonic oscillator frequency 
is selected. In our work we have circumvented these dif- 
ficulties by closely following the biquadratic potential to 
determine the corresponding parameters of the Gaussian 
model, as noted in Sec. |V A| A fixed harmonic oscilla- 
tor frequency has been used throughout the calculation 
for a particular charge stability diagram, which corre- 
sponds to a point in the (1,1) regime where the potential 
wells are symmetric. The motivation of this constraint is 
that the biquadratic model potential, with loq fixed, has 
been shown to well describe the experiment with much 
less n umber of independent variables than the Gaussian 
modelP^im We believe that in order to gain physical in- 
sight it is sufficient to manipulate the few most important 
variables. A precise treatment for the Gaussian model, 
with the harmonic oscillator frequencies appropriately 
vary with the deformation of the potential well, gives 
results that are not quite different from the calculations 
keeping loq fixed in comparison to experiments. How- 



ever, it is useful to briefly study how the charge stability 
diagram changes as the harmonic oscillator frequency is 
varied. 

FigJT^a) shows the identical result as that shown in 
Fig. fTofa). The harmonic oscillator frequency (energy) 
has been chosen at the locus denoted by a white cross: 
hujQG — 6.487 meV. In the calculation of Fig. [l2^b), 
we artificially multiply the harmonic oscillator energy by 
a factor of 1.5, thus further localizing the Fock-Darwin 
wave functions but leaving the confinement potential un- 
changed. Beside an overall shift in energy (due to the 
larger Hujq), the (1,1) component of the charge stabil- 
ity diagram clearly expands its area, meaning that in 
this case the on-site Coulomb interaction is increased. 
In Fig. [l2^c) , the charge stability diagram is calculated 
using Fock-Darwin states whose harmonic oscillator fre- 
quencies are carefully calculated according to the precise 
form of Gaussian model at each given point on the charge 
stability diagram. We see that the lower left components 
(N < 2) are quite similar to that of Fig. 12 'a), while 
for N > 2 the spatial variation of the on-site Coulomb 
interaction U is evident. The shape of the (1,1) com- 
ponent of the charge stability diagram is abnormal: it 
shrinks its area in proximity to the N = 3 block, while its 
boundaries with the N = 1 block remain approximately 
the same. As a consequence, pairs of the boundaries 
of the (1,1) regime are not parallel, in contrast to what 
observed in other calculations and in most experiments. 
This means that the case considered here, varying the 
harmonic oscillator frequency of the Fock-Darwin states 
as the potential wells are changing, is less relevant to ex- 
periments that we have considered. The good correspon- 
dence between Fig. 12 'a) and (c) means that the approxi- 
mation that we have employed is reasonable; Rather, the 
hoj artificially chosen for Fig. (l2fcb) is too large to accu- 
rately reflect the actual solution. We therefore conclude, 
since the present calculations already have quantitative 
correspondence with experiments, that the method we 
have used in previous sections using Fock-Darwin states 
with fixed harmonic oscillator frequency is most suitable 
for our work. 

In Fig. [i"27d) we show a calculation with parameters 
same as Fig. [12fa) except that there is a finite external 
magnetic field B = 6 T. The reason that we group this 
result with others discussed in this section is that the 
magnetic field can be regarded to have very similar ef- 
fect as increasing the harmonic oscillator frequency, as 
can be seen by similarities between Fig. 12 ^b) and (d). 
In fact, this can be straightforwardly understood from 



Eq. (12 1: A finite magnetic field leads to a finite Larmor 
frequency ui c , which consequently enters the energy of 
the Fock-Darwin states in a similar way as the harmonic 
oscillator frequency (at dot j) uij. Therefore, a magnetic 
field effectively localizes the electron wave function, lead- 
ing to stronger on-site Coulomb interactions and weaker 
inter-site Coulomb interactions. 
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FIG. 13: Illustration of the connections between different 
components discussed in this paper as well as Refs. [50] and 

ED 



VI. CONCLUSION 

In this paper, along with our previous publications 
Refs. [50] and [51] we have developed a theoretical frame- 
work treating the coupled double quantum dot system 
with a generalized form of Hubbard model. The aim of 
our work is to establish conceptual and quantitative con- 
nections between experiments, existing theoretical ap- 
proaches and the new theory (see Fig. 13), thus present- 



ing a unified and transparent description of our current 
understanding of the double quantum dot system and 
allowing possibilities for studying more complicated sit- 
uations in future work. The central object of these stud- 
ies, which is invariably present in all these connections, 
is the charge stability diagram, the visualization of elec- 
tron configurations in the double dot system. Given the 
essential role of semiconductor double quantum systems 
in solid state quantum computer architectures, our use 
of the generalized Hubbard model in describing the sys- 
tem is natural since the necessary entanglement effects 
are embedded in the interaction terms of the Hubbard 
model. 

Experimentally, the charge stability diagram can be 
directly measured by a quantum point contact,^ and it 
is therefore a genuine representation of electron states 
in the quantum dot system. It contains two kinds of 
information: the classical Coulomb effect and quantum 
fluctuations. The classical aspect of the change stability 
diagram is well understood using the capacitance circuit 
modelPHU] On the other hand, microscopic calculations 
are usuall y carried out for other aspects of the quantum 
dot systeirPsESmHlZ] onr y rarely for the charge sta- 
bility diagram.^ This means that our understanding of 
the quantum aspect of the charge stability diagram has 
largel y been limited to the quantum-mech anical tw o-level 
modelJSli^riHl w ith only a few except ions EMUSil 



The generalized Hubbard model (see Sec. Ill), intro- 
duced in Ref. 1501 has been naturally written down with 
only the simplest symmetry considerations in mind. It 
includes all terms allowed by symmetry, which therefore 
means that all possible quantum fluctuation effects are 
included. On one hand, when the quantum fluctuations 



are completely suppressed, the model is mapped exactly 
to the capacitance circuit model. On the other hand, 
quantum fluctuations perturb the charge s tabili ty dia- 
gram, which is observed both in experiment j 32 * 33 ^ and in 
theoretical calculations (see Sec. IV). The intrinsic as- 
sumption of the capacitance circuit model is that the 
electrostatic energy of a given system are functions of in- 
teger electron occupancies on the dots. We have argued 
that this assumption eliminates the possibility of directly 
including quantum fluctuations in the capacitance cir- 
cuit model, as the electron occupancies cease to be good 
quantum numbers in the presence of quantum fluctua- 
tions. Therefore, the generalized Hubbard model is the 
most direct quantum generalization of the capacitance 
circuit model. 

The generalized Hubbard model contains several pa- 
rameters. Although their values are not determined by 
symmetry considerations, they can not be arbitrarily cho- 
sen, otherwise unphysical charge stability diagrams may 
occur (see Appendix [B]) . On one hand, we calculate all 
parameters of the generalized Hubbard model through 
the microscopic theory, detailed in Sec. [H] The micro- 
scopic calculation involves a numerical solution to the 
many-electron Schrodinger equation under some model 
confinement potential. The electron wave functions are 
approximated by linear combinations of the Fock-Darwin 
states through the configuration interaction method. On 
the other hand, some of the parameters can be extracted 
from experimentally measured charge stability diagrams, 
provided that the resolution is sufficiently high. There- 
fore the generalized Hubbard model acts as a bridge 
which quantitatively connects the experiments and the 
microscopic theory. These considerations, together with 
the connections to the capacitance model, are shown 
symbolically in Fig. |13| Fig. [13] implies that our work in 
this paper as well as Refs. [511] and [ST] present a coherent 
and complete theoretical picture of our current under- 
standing of the double quantum dot system including all 
effects of quantum fluctuations and entanglement. 

In our previous works JsSlSU we have introduced the gen- 
eralized Hubbard model, elaborated the motivation and 
various simplifications, and applied it to explain several 
specific experiments. We consider the present paper as 
the concluding one of this three-paper series studying the 
charge stability diagram of double quantum dot system. 
In other words, the goal of this paper is to complete the 
theoretical framework that we have developed, providing 
detailed in forma tion that has not been covered by previ- 
ous papers J^lliU thus forming a comprehensive and coher- 
ent picture together with the previous works. First, in 
Refs. l50lST1 we have focused on one of the most important 
quantum parameters in the generalized Hubbard model, 
namely the tunnel coupling, but other quantum param- 
eters are important as well. In this paper, we have stud- 
ied the effect of all quantum parameters on the charge 
stability diagram (Sec. IV), and several cases with ex- 
treme parameters are shown (Appendix [B]) . Second, in 
Ref. [50] we have only used the biquadratic potential as 
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the model confinement potential in the microscopic cal- 
culations. Although in Ref . ISTl results from both the bi- 
quadratic and Gaussian models are shown, a direct com- 
parison between the two models is missing. In this paper, 
we have directly compared the results of the charge sta- 
bility diagram calculated using these two model poten- 
tials (Sec. VA). We have found that since the Gaussian 



model has narrower potential wells, the electron states 
are more localized, leading to a stronger on-site Coulomb 
interaction and weaker inter-site quantum fluctuations. 
Third, in all previous works we have retained the lowest- 
lying orbital only. In order to justify this simplification, 
in this paper we have studied the effect of higher-lying or- 
bitals on the one -orbi tal subspace that we have previously 
considered (Sec. VB). We have found that the additional 



p-orbitals lead to a small renormalization of Coulomb in- 
teraction parameters while leaving the main physics un- 
changed. This justifies our practice of focusing on the 
one-orbital subspace. Last, we have studied the role of 
the harmonic oscillator frequency in the Gaussian model 
and the effect of external magnetic field (Sec. VCl. We 
have found that the electron wave functions with fixed 
value of harmonic oscillator frequency is most suitable 
for comparison with experiments. The influence of an 
external magnetic field can be understood as effectively 
increasing the harmonic oscillator frequency in the Fock- 
Darwin wave function. Both of these are not covered by 
our previous works. 

As discussed in Sec. [TlJ we have neglected the inter- 
action of the quantum dot spin qubit with the envi- 
ronment because environmental decoherence, which is 
a huge subject by itself, is well beyond the scope of 
the current paper dealing with the qubit control and 
charge stability diagram in semiconductor quantum dots. 
Various environment factors contribute to the d ecoher- 
ence process, which include the charge noise J2U221 arising 
from the electromagnetic fluctuations in the surround- 
ing gates and th e dynamics in unintentional background 
impurities) 64 ^ 65 ! the hyperfine cou pling to the surround- 
ing nuclear spin bathj 4 ^ 4 - 2 1 43 1 62 1 631 spectral diffusion of the 
electron spin due to slow nuclear spin dynamics arising 
from the nuclear spin flip-fl ops due to dipolar coupling 
within the nuclear spin bathp2MMl th e coupling of elec- 
tron spins to unintentional background inuyurity spins in 
the environment P^and phonon modesPSlE2l In many sit- 
uations, the most important environmental factor is the 
effect of unintentionally introduced background charged 
impurities, since depending on the proximity between 
the impurities and the qubit, the interaction with the 
impurities could be of the same order of magnitude as 
the on-dot electronic interactions as both are mediated 
by the long-range Coulomb couplingPUl The environ- 
mental nuclear spin bath affects the electron spin both 
through a "direct" hyperfine interaction between the spin 
qubit and the nuclei, and an "indirect" coupling due to 
the spectral diffusion process resulting from the dipo- 
lar intera ction between the nuclei. Recent experimental 
progresse d 18 ^ 24 ! have demonstrated control over the nu- 



clear spin bath, making the nuclear spin bath less de- 
structive and the effect of the impurities more prominent. 
In some systems such as Si or C, the direct hyperfine cou- 
pling can be eliminated by using isotropic purification of 
the nuclear elements! py Phonon modes are the least im- 
portant environmental factor for spin qubits working at 
low (~ 100 mK) temperatures: Firstly, the coupling be- 
tween the spin qubit and the phonon modes are orders of 
magnitude weaker than the effect of charged impurities 
and nuclear spin bath mentioned above. Secondly, the 
effect of phonon modes can be suppressed by cooling the 
system down to very low temperatures, while the effects 
due to charged impurities and the nuclear spin bath re- 
main important even at low temperatures. In this paper, 
we focus on the electronic interaction within the quan- 
tum dot system and study its consequence on the charge 
stability diagram, and quantum decoherence due to en- 
vironmental factors is beyond the scope of the current 
work. 

Our work suggests several possible directions for future 
studies. First, it is useful to go beyond the model poten- 
tial and Hund-Mulliken approximation, and develop a 
first-principles Poission-Schrodinger approach for a real- 
istic calculation of quantum dot confinement. This in- 
volves both extracting the detailed form of confinement 
potentials from experimentally fabricated devices, and 
an exact treatment of Schrodinger equation under this 
realistic confinement potential. However, this would re- 
quire a precise experimental control of disorder and im- 
purities during the fabrication, as well as the knowledge 
of the underlying Fermi sea. Moreover, the calculation 
expenses are much higher than what we have done in 
this work, although they are not necessarily prohibitive 
since modern numerical techniques, e.g. density func- 
tional theory, have already demonstrated their power in 
successfully explaining many aspects of complex mate- 
rials as well as nanoscale systems. Second, although we 
have argued that the additional higher-lying orbitals pose 
no more than a renormalization of Coulomb parameters 
in the single-orbital subspace, the issue for the full multi- 
orbital mult i- electron problem may not yet be settled. A 
detailed study of the realistic multi-orbital multi-electron 
problem would in particular greatly benefit our under- 
standing of the silicon quantum dot systerrPH as this 
multi-electron effect may be one of the obstacles of real- 
izing and manipulating the singlet/triplet qubit in silicon 
double quantum dot system. Last, we have focused on 
equilibrium properties of an isolated quantum dot system 
where the total electron number and the total spin are 
conserved. However, non-equilibrium effects are inter- 
twined with the measurement of the charge stability di- 
agram as well as the manipulation of qubits. The charge 
stability diagram is measured through a tunnelling cur- 
rent, which itself is a non-equilibrium property. In Ref. 1811 
a study of the charge stability diagram in terms of the 
microscopic conductance(admittance) is given, and much 
of our theory can also be applied along this direction. Of 
course, such a dynamical theory would be much more 
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complicated than what we have studied here, since the 
excited states must be considered in all non-equilibrium 
calculations. Moreover, the coupling of the quantum dot 
system to its environment, which has been neglected in 
this paper, is by no means insignificant. Therefore, the 
assumptions that we have made in this paper as well as 
Refs. 50 51, that the system is isolated and in equilib- 
rium, should be lifted in future studies by extending the 
generalized Hubbard model. The theoretical approach 
that we have presented here provides a foundation for 
these future research directions. 
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Appendix A: Matrix form of the generalized 
Hubbard model 

The full Hamiltonian is a 16 x 16 matrix. Because 
the total particle number (N) and total spin (S z ) are 



conserved, it appears in a block-diagonal form, as follows. 

(1) N = 0, S z = 0. On the basis |0), 

H x = 0. (A-l) 

(2) N = 1, S z = §. On the basis |c} t |0),4 t |0)}, 

M"" 1 -^" -„+«,)■ (A - 2) 

(3) N = 1, S z = -i. On the basis |c f u |0) |0)|, 

":<={-'"-< Eb -»;' Eb )- 

(4) N = 2, S z = 1. On the basis cLcL |0) (labeled as 

Hi = -m - fi 2 + U 12 - J e + 2E B . ( A-4) 

(5) N = 2, S z = 0. On the basis 

{AAi i°) . 4 t c u i°) . AAi i°) > 4t c l 1°)}' 



-^1-^2+^12 Je -(t + Jtl) -{t+Jtt) 

j, _ \ J e -A*i- M2 + 012 -(t + Jti) -(t + Jti) I , \ - , 

W5 ~' -(t + Jti) -(t + J a ) -2/ii + f/a J p I {A ~ b) 



-(t + J t2 ) -(t + J t2 ) J p -2fi 2 + U: 



2 



Note that the linear combinations (4f c 2,|. — 4t c u) (labeled as T ) and f4f c 24. + 4t c u) (labeled as 5) 
decouple the To state with energy — /ii — fi 2 + f 12 — Je, thus forming a triplet which degenerate at zero magnetic field. 
The remaining 3x3 matrix reads 

/ -/il - M2 + 012 + Je -V2(t + J tl ) -\/2(t + J t2 ) \ 

#5 = -V2(t+J tl ) -2 Ml + f/i J p . (A-6) 

\ -V2(t+J t2 ) J P -2/z 2 + 02 J 

(6) N = 2, S z = -1. On the basis cj.ct |0) (labeled as T_), 

#6 = -Mi - M2 + 012 - Je - 2E B . ( A-7) 



(7) iV = 3, 5* = §. On the basis {4-|-4r c u |0) ,c\ x c\ x c\i I )}' 



^ _ I -2m - H2 + Ul + 2012 - Je + #B -(i+Jil+J t2 ) ] / a 

' I _(t + j tl+ j a ) - 2f i 2 -fi 1 + U 2 +2U 12 -J e +E B l- 



(8) JV = 3, S z = -§. On the basis |4t c u4i 1°) > 4t c H c L I )}' 



„ _ , -2ah - M2 + 01 + 20i2 - J e - + + Je) 1 , \ ,,, 

' N 1 -(t+ Ju + J* 2 ) -2^2 - Mi + 02 + 2C/12 -J e -E F 1 ' 



(9) AT = 4, (5 Z = 0. On the basis 4t c 2t c Ii c 24. |0), Appendix B: Several extreme cases of the 



H g = -2/n - 2^ 2 + U x + U 2 + 4U 12 - 2 J e . (A-10) 



parameters 

In this section we present several extreme cases of pa- 
rameters. As discussed above, J e < Ui 2 < U in a typical 
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FIG. 14: Charge stability diagram calculated at U = 6.1 meV, 
Urn = 7 meV, t = J e = J p = J t = 0. 



physical system. However it remains interesting to study 
the cases with U X2 > U or J e > t/i2- First, for complete- 
ness of the understanding of the charge stability diagram, 
these parameter regimes need be explored. Indeed, as 
shall be seen below, the topology of the charge stability 
diagram changes. Second, although the systems we are 
currently studying impose special limitations of param- 
eters, this in general need not be true. For example, in 
our study J e <C t. However, in a lattice problem the 



orthogonality of atomic orbitals means that t = while 
J e is still finite. Third, interesting physical phenomena 
happen in some of the parameter regimes, such as the 
U12 > U case discussed in Ref. which is at least of 
theoretical interest. 

Fig. [14] shows a typical case for U X2 > U. The coor- 
dinates (/i2, Hi) of the triple points are given by A(0, 0), 
B(U, U), B'(2U 12 , 2U 12 ), A'{U + 2U X2 , U + 2U X2 ). The 
(1,1) component is completely ruled out from the charge 
stability diagram since the system always favors either 
(0, 2) or (2, 0) state which have smaller Coulomb repul- 
sion than the (1, 1) state. This charge stability diagram is 
not observed experimentally, which means that it is rea- 
sonable to set U12 < U for currently implemented double 
quantum dot systems. 

In Fig. [l5] we show the case with J e > Ui 2 (but still 
J e < U). Fig. 15 a) shows the J e = U± 2 result, which can 
be understood as the limit of Fig. [4] Fig. [Ts^b) shows 
the result of Ui 2 < J e < U. The topology fundamen- 
tally changes in that the (1,1) triplet state now shares a 
phase boundary directly with the (0,0) state. The coor- 
dinates (^*2>Mi) °f the triple points are A(Ui 2 — J e) 0), 
B(0,Ui2 — J e ), and C(U,Ui 2 — J e ), with others found 
by symmetry. Interestingly the length of line segment 
AB is still V2(J e — U12) (which can be expressed as 
\/2 \ Ui2 — J e \ formally similar to that of Fig. GJ. 



1 R. P. Feynman, Int. J. Theor. Phys. 21, 467 (1982); Found. 
Phys. 16, 507 (1986); D. Deutsch, Proc. R. Soc. London, 
Ser. A 400, 97 (1985) 

2 P. W. Shor, in Proceedings of the 35th Annual Symposium 
on the Foundations of Computer Science (IEEE Press, Los 
Alamitos, CA, 1994), pp. 124133; A. Ekert and R. Jozsa, 
Rev. Mod. Phys. 68, 733 (1996); L. K. Grover, Phys. Rev. 
Lett. 79, 325 (1997); A. Steane, Rep. Prog. Phys. 61, 117 
(1998); C. H. Bennett and D. P. DiVincenzo, Nature (Lon- 
don) 404, 247 (2000). 

3 L. M. K. Vandersypen, M. Steffen, G. Breyta, C. S. Yan- 
noni, R. Cleve, and I. L. Chuang, Phys. Rev. Lett. 85, 5452 
(2000); L. M. K. Vandersypen, M. Steffen, G. Breyta, C. 
S. Yannonil, M. H. Sherwood, and I. L. Chuang, Nature 
(London) 414, 883 (2001) 

4 T. Sleator and H. Weinfurter, Phys. Rev. Lett. 74, 4087 
(1995); Q. A. Turchette, C. J. Hood, W. Lange, H. 
Mabuchi, and H.J. Kimble, ibid. 75, 4710 (1995). 

5 J. I. Cirac and P. Zoller, Phys. Rev. Lett. 74, 4091 (1995); 

C. Monroe, D. M. Meekhof, B. E. King, W. M. Itano, and 

D. J. Wineland, ibid. 75, 4714 (1995); C. A. Sackett, D. 
Kielpinski, B. E. King, C. Langer, V. Meyer, C. J. Myatt, 
M. Rowe, Q. A. Turchette, W. M. Itano, D. J. Wineland, 
and C. Monroe, Nature (London) 404, 256 (2000). 

6 D. Loss and D. P. DiVincenzo, Phys. Rev. A 57, 120 
(1998). 

7 D. P. DiVincenzo, D. Bacon, J. Kempe, G. Burkard, and 
K. B. Whaley, Nature 408, 339 (2000). 

8 J. Levy, Phys. Rev. Lett. 89, 147902 (2002). 



9 F. H. L. Koppens, C. Buizert, K. J. Tielrooij, I. T. Vink, 
K. C. Nowack, T. Meunier, L. P. Kouwenhoven, and L. M. 
K. Vandersypen, Nature 442, 766 (2006). 

B. E. Kane, Nature 393, 133 (1998). 

1 R. Vrijen, E. Yablonovitch, K. Wang, H. W. Jiang, A. 
Balandin, V. Roychowdhury, T. Mor, and D. DiVincenzo, 
Phys. Rev. A 62, 012306 (2000). 

2 M. Friesen, P. Rugheimer, D. E. Savage, M. G. Lagally, D. 
W. van der Weide, R. Joynt, and M. A. Eriksson, Phys. 
Rev. B 67, 121301(R) (2003). 

3 J. M. Taylor, H.-A. Engel, W. Diir, A. Yacoby, C. M. Mar- 
cus, P. Zoller, and M. D. Lukin, Nat. Phys. 1, 177 (2005). 

4 J. M. Kikkawa and D. D. Awschalom, Phys. Rev. Lett. 80, 
4313 (1998). 

5 S. Amasha, K. MacLean, Iuliana P. Radu, D. M. Zumbuhl, 
M. A. Kastner, M. P. Hanson, and A. C. Gossard, Phys. 
Rev. Lett. 100, 046803 (2008). 

6 F. H. L. Koppens, K. C. Nowack, and L. M. K. Vander- 
sypen, Phys. Rev. Lett. 100, 236802 (2008). 

7 C. Barthel, J. Medford, C. M. Marcus, M. P. Hanson, and 
A. C. Gossard, Phys. Rev. Lett. 105, 266808 (2010). 

8 H. Bluhm, S. Foletti, I. Neder, M. Rudner, D. Mahalu, V. 
Umansky, and A. Yacoby, Nat. Phys. 7, 109 (2011). 

9 A. M. Tyryshkin, S. A. Lyon, A. V. Astashkin, and A. M. 
Raitsimring, Phys. Rev. B 68, 193207 (2003). 

10 A. Morello, J. J. Pla, F. A. Zwanenburg, K. W. Chan, K. 
Y. Tan, H. Huebl, M. Mottonen, C. D. Nugroho, C. Yang, 
J. A. van Donkelaar, A. D. C. Alves, D. N. Jamieson, C. 
C. Escott, L. C. L. Hollenberg, R. G. Clark, and A. S. 



18 



10 



> 

s 5 



-5 



'■_ (2,0) 


1 1 


i i 1 i i i 
(2,1) 


■ J • ■ ■ _ 
(2,2) ~ 












: (1,0) 


(1,1)7 




(1,2) J 










: (o,o) 

i i i i 


(0,1) 
i i i i I 


i i 


(0,2) : 
i I i i i " 



(a) 



-5 5 10 

ju 2 (meV) 



10 



I 1 1 1 1 I 
- (2,0) 



1 1 1 I 1 
(2,1) 



(2,2) 



(1,0) 



A 



(1,2) : 



(0,0) 



(0,1) 



(0,2) 



(b) 



-5 5 10 

ju 2 (meV) 



FIG. 15: Charge stability diagram calculated at U = 6.1 meV, 
U12 = 2.5 meV, (a) J e = 2.5 meV = U 12 (b) J e = 4 meV > 
Ui 2 . t = J p = J t = 0. 



Dzurak, Nature (London) 467, 687 (2010). 

21 C. B. Simmons, J. R. Prance, B. J. Van Bael, T. S. Koh, 
Z. Shi, D. E. Savage, M. G. Lagally, R. Joynt, M. Friesen, 
S. N. Coppersmith, and M. A. Eriksson, Phys. Rev. Lett. 
106, 156804 (2011). 

22 M. Xiao, M. G. House, and H. W. Jiang, Phys. Rev. Lett. 
104, 096801 (2010). 

23 J. R. Petta, A. C. Johnson, J. M. Taylor, E. A. Laird, A. 
Yacoby, M. D. Lukin, C. M. Marcus, M. P. Hanson, A. C. 
Gossard, Science 309, 2180 (2005). 

24 S. Foletti, H. Bluhm, D. Mahalu, V. Umansky, and A. 
Yacoby, Nat. Phys. 5, 903 (2009). 

25 E. A. Laird, J. M. Taylor, D. P. DiVincenzo, C. M. Marcus, 
M. P. Hanson, and A. C. Gossard, Phys. Rev. B 82, 075403 
(2010). 

26 I. van Weperen, B. D. Armstrong, E. A. Laird, J. Medford, 
C. M. Marcus, M. P. Hanson, and A. C. Gossard, Phys. 
Rev. Lett. 107, 030506 (2011). 

27 E. P. Nordberg, G. A. Ten Eyck, H. L. Stalford, R. P. 
Muller, R. W. Young, K. Eng, L. A. Tracy, K. D. Childs, 
J. R. Wendt, R. K. Grubbs, J. Stevens, M. P. Lilly, M. 
A. Eriksson, and M. S. Carroll, Phys. Rev. B 80, 115331 
(2009). 

28 S. Das Sarma, R. de Sousa, X. Hu, and B. Koiller, Sol. St. 
Comm. 133, 737 (2005). 

29 D. Culcer, L. Cywihski, Q. Li, X. Hu, and S. Das Sarma, 



Phys. Rev. B 80, 205302 (2009); 82, 155312 (2010). 

30 S. Goswami, K. A. Slinker, M. Friesen, L. M. McGuire, J. 
L. Truitt, C. Tahan, L. J. Klein, J. O. Chu, P. M. Mooney, 
D. W. van der Weide, R. Joynt, S. N. Coppersmith, and 
M. A. Eriksson, Nat. Phys. 3, 41 (2007). 

31 M. G. Borselli, R. S. Ross, A. A. Kiselev, E. T. Croke, K. 
S. Holabird, P. W. Deelman, L. D. Warren, I. Alvarado- 
Rodriguez, I. Milosavljevic, F. C. Ku, W. S. Wong, A. E. 
Schmitz, M. Sokolich, M. F. Gyure, and A. T. Hunter, 
Appl. Phys. Lett. 98, 123118 (2011). 

32 C. B. Simmons, M. Thalakulam, B. M. Rosemeyer, B. J. 
Van Bael, E. K. Sackmann, D. E. Savage, M. G. Lagally, R. 
Joynt, M. Friesen, S. N. Coppersmith, and M. A. Eriksson, 
Nano Lett. 9 3234 (2009). 

33 N. S. Lai, W. H. Lim, C. H. Yang, F. A. Zwanenburg, 
W. A. Coish, F. Qassemi, A. Morello, and A. S. Dzurak, 
e-print |arXiv:1012.1410| (unpublished) . 

34 C. Livermore, C. H. Crouch, R. M. Westervelt, K. L. 
Campman, A. C. Gossard, Science 274, 1332 (1996). 

35 J. M. Elzerman, R. Hanson, J. S. Greidanus, L. H. Willems 
van Beveren, S. De Franceschi, L. M. K. Vandersypen, 
S. Tarucha, and L. P. Kouwenhoven, Phys. Rev. B 67, 
161308(R) (2003). 

36 F. Hofmann, T. Heinzel, D.A. Wharam, J.P. Kotthaus, G. 
Bohm, W. Klein, G. Trankle, and G. Weimann, Phys. Rev. 
B 51, 13872 (1995). 

37 L. DiCarlo, H. J. Lynch, A. C. Johnson, L. I. Childress, K. 
Crockett, C. M. Marcus, M. P. Hanson, and A. C. Gossard, 
Phys. Rev. Lett. 92, 226801 (2004). 

38 V. W. Scarola, K. Park, and S. Das Sarma, Phys. Rev. 
Lett. 93, 120503 (2004); V. W. Scarola and S. Das Sarma, 
Phys. Rev. A 71, 032340 (2005). 

39 Q. Li, L. Cywihski, D. Culcer, X. Hu, and S. Das Sarma, 
Phys. Rev. B 81, 085313 (2010). 

40 J. M. Taylor, J. R. Petta, A. C. Johnson, A. Yacoby, C. M. 
Marcus, and M. D. Lukin, Phys. Rev. B 76, 035315 (2007) 

41 J. R. Petta, J. M. Taylor, A. C. Johnson, A. Yacoby, M. D. 
Lukin, C. M. Marcus, M. P. Hanson, and A. C. Gossard, 
Phys. Rev. Lett. 100, 067601 (2008). 

42 D. J. Reilly, J. M. Taylor, J. R. Petta, C. M. Marcus, M. 
P. Hanson, and A. C. Gossard, Science 321, 817 (2008). 

43 D. J. Reilly, J. M. Taylor, J. R. Petta, C. M. Marcus, M. P. 
Hanson, and A. C. Gossard, Phys. Rev. Lett. 104, 236802 
(2010). 

44 C. Barthel, D. J. Reilly, C. M. Marcus, M. P. Hanson, and 
A. C. Gossard, Phys. Rev. Lett. 103, 160503 (2009). 

45 C. Barthel, M. Kjaergaard, J. Medford, M. Stopa, C. M. 
Marcus, M. P. Hanson, and A. C. Gossard, Phys. Rev. B 
81 161308(R), (2010). 

46 C. W. J. Beenakker, Phys. Rev. B 44, 1646 (1991) 

47 W. G. van der Wiel, S. De Franceschi, J. M. Elzerman, 
T. Fujisawa, S. Tarucha, L. P. Kouwenhoven, Rev. Mod. 
Phys. 75, 1 (2003). 

48 R. Hanson, L. P. Kouwenhoven, J. R. Petta, S. Tarucha, 
L. M. K. Vandersypen, Rev. Mod. Phys. 79, 1217 (2007). 

49 D. Schroer, A. D. Greentree, L. Gaudreau, K. Eberl, L. C. 
L. Hollenberg, J. P. Kotthaus, and S. Ludwig, Phys. Rev. 
B 76, 075306 (2007). 

50 S. Yang, X. Wang, and S. Das Sarma, Phys. Rev. B 83, 
161301(R) (2011). 

51 S. Das Sarma, X. Wang, and S. Yang, Phys. Rev. B 83, 
235314 (2011). 

52 J. R. Petta, A. C. Johnson, C. M. Marcus, M. P. Hanson, 
and A. C. Gossard, Phys. Rev. Lett. 93, 186802 (2004) 



19 



53 T. Hatano, M. Stopa, and S. Tarucha, Science 309, 268 
(2005). 

54 A. K. Hiittel, S. Ludwig, H. Lorenz, K. Eberl, and J. P. 
Kotthaus, Phys. Rev. B 72, 081310 (2005). 

55 M. Pioro-Ladriere, M. R. Abolfath, P. Zawadzki, J. La- 
pointe, S. A. Studenikin, A. S. Sachrajda, and P. Hawry- 
lak, Phys. Rev. B 72, 125307 (2005). 

56 L.-X. Zhang, D. V. Melnikov, and J.-P. Leburton, Phys. 
Rev. B 74, 205306 (2006). 

57 C. A. Stafford and S. Das Sarma, Phys. Rev. Lett. 72, 
3590 (1994). 

58 R. Kotlyar, C. A. Stafford, and S. Das Sarma, Phys. Rev. 
B 58, R1746 (1998); 58, 3989 (1998); C. A. Stafford, R. 
Kotlyar, and S. Das Sarma, ibid. 58, 7091 (1998). 

59 J. H. Jefferson and W. Hausler, Phys. Rev. B 54, 4936 
(1996). 

60 L. Gaudreau, S. A. Studenikin, A. S. Sachrajda, P. Za- 
wadzki, A. Kam, J. Lapointe, M. Korkusinski, and P. 
Hawrylak, Phys. Rev. Lett. 97, 036807 (2006). 

61 M. Korkusinski, I. P. Gimenez, P. Hawrylak, L. Gaudreau, 
S. A. Studenikin, and A. S. Sachrajda, Phys. Rev. B 75, 
115301 (2007). 

62 W. M. Witzel and S. Das Sarma, Phys. Rev. B 74, 035322 
(2006) 

63 L. Cywifiski, W. M. Witzel, and S. Das Sarma, Phys. Rev. 
Lett. 102, 057601 (2009); Phys. Rev. B 79, 245314 (2009) 

64 I. P. Gimenez, C.-Y. Hsieh, M. Korkusinski, and P. Hawry- 
lak, Phys. Rev. B 79, 205311 (2009) 

65 N. T. T. Nguyen and S. Das Sarma, Phys. Rev. B 83, 
235322 (2011) 

66 M. J. Storcz, U. Hartmann, S. Kohler, and F. K. Wilhelm, 
Phys. Rev. B 72, 235321 (2005) 

67 V. N. Stavrou and X. Hu, Phys. Rev. B 72, 075362 (2005) 

68 P. Stano and J. Fabian, Phys. Rev. Lett. 96, 186602 (2006) 

69 D. Harbusch, D. Taubert, H. P. Tranitz, W. Wegscheider, 
and S. Ludwig, Phys. Rev. Lett. 104, 196801 (2010) 

70 X. Hu, Phys. Rev. B 83, 165322 (2011) 

71 E. Artacho and L. Milans del Bosch, Phys. Rev. A 43, 5770 
(1991). 

72 X. Hu and S. Das Sarma, Phys. Rev. A 61, 062301 (2000). 

73 J. C. Slater, Phys. Rev. 49, 537 (1936). 

74 P. W. Anderson, Phys. Rev. 115, 2 (1959). 

75 J. Kanamori, J. Phys. Chem. Solids 10, 87 (1959); Prog. 



Theor. Phys. 30, 275 (1963). 

76 J. E. Hirsch, Physica C 199, 305 (1992); J. E. Hirsch and 
F. Marsiglio, Phys. Rev. B 62, 15131 (2000). 

77 J. Hubbard, Proc. R. Soc. London, Ser. A 276, 238 (1963); 
277, 237 (1964); 281, 401 (1964). 

78 J. Zaanen and M. Oles, Ann. Phys. (Leipzig) 5, 224 (1996) 

79 M. Imada, A. Fujimori, Y. Tokura, Rev. Mod. Phys. 70, 
1039 (1998). 

80 A. Liebsch, Phys. Rev. Lett. 91, 226401 (2003); A. Koga, 
N. Kawakami, T. M. Rice, and M. Sigrist, Phys. Rev. Lett. 
92, 216402 (2004); P. Werner and A. J. Millis, Phys. Rev. 
B 74, 155107 (2006) 

81 A. Cottet, C. Mora, and T. Kontos, Phys. Rev. B 83, 
121311(R) (2011). 

82 G. Burkard, D. Loss, and D. P. DiVincenzo, Phys. Rev. B 
59, 2070 (1999). 

83 X. Hu and S. Das Sarma, Phys. Rev. A 64, 042312 (2001). 

84 R. de Sousa, X. Hu, and S. Das Sarma, Phys. Rev. A 64, 
042307 (2001). 

85 I. Puerto Gimenez, M. Korkusinski, and P. Hawrylak, 
Phys. Rev. B 76, 075336 (2007). 

86 E. Nielsen, R. W. Young, R. P. Muller, and M. S. Carroll, 
Phys. Rev. B 82, 075319 (2010) 

87 E. Nielsen and R. P. Muller, e-print |arXiv:1006.2735l (un- 
published) 

88 M. Helle, A. Harju, and R. M. Nieminen, Phys. Rev. B 72, 
205329 (2005) 

89 J. Pedersen, C. Flindt, N. A. Mortensen, and A.-P. Jauho, 
Phys. Rev. B 76, 125323 (2007) 

90 J. G. Pedersen, C. Flindt, A.-P. Jauho, and N. A. 
Mortensen, Phys. Rev. B 81, 193406 (2010). 

91 X. Hu and S. Das Sarma, Phys. Rev. Lett. 96, 100501 
(2006) 

92 M. Stopa, and C. M. Marcus, Nano Lett., 8, 1778 (2008). 

93 W. M. Witzel, R. de Sousa, and S. Das Sarma, Phys. Rev. 
B 72, 161306(R) (2005) 

94 W. M. Witzel and S. Das Sarma, Phys. Rev. B 77, 165319 
(2008) 

95 W. M. Witzel, M. S. Carroll, A. Morello, L. Cywinski, and 
S. Das Sarma, Phys. Rev. Lett. 105, 187602 (2010) 

96 X. Wang and A. J. Millis, Phys. Rev. B 81, 045106 (2010). 



